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Geomorphometry  and  Quality  Control 

Peter  L.  Guth 

Department  of  Oceanography,  U.S.  Naval  Academy, 

Annapolis,  Maryland  21402-5026,  U.S. A. 
e-mail:  pguth  @  usna.  edit 

Abstract.  Digital  elevation  models  (DEMs)  cover  a  wide  range  of  scales,  and  allow 
statisticalanalysis  of  geomorphometric  parameters.  At  global  or  continental  scale,  DEMs 
covering  rectangular  quadrangles  can  be  considered  random  samples.  Average  quadrangle 
values  can  be  compared  at  different  DEM  scales,  and  for  different  physiographic 
provinces.  Three  independent  variables  provide  valuable  descriptions  of  terrain:  average 
elevation,  average  slope,  and  the  degree  of  terrain  organization.  DEMs  with  spacings  of 
30"  (global  coverage),  3"  (continental  United  States),  and  30  m  (local  United  States) 
provide  almost  perfect  correlations  for  average  quadrangle  elevation  and  slope,  although 
the  slope  values  increase  as  the  DEM  spacing  decreases.  The  degree  of  terrain  organization 
also  correlates  across  DEM  scales,  but  with  lower  correlation  coefficients,  especially  in 
areas  of  lower  relief.  Terrain  variables  computed  from  10  m  and  30  m  USGS  Level  2 
DEMs  are  essentially  identical.  Slope  algorithms  perform  differently  in  different 
physiographic  provinces,  and  the  aspect  algorithm  performs  poorly  in  low  relief  areas. 
Geomorphometric  analysis  can  provide  a  rapid  and  effective  assessment  of  DEM  quality 
control,  and  should  be  integrated  into  the  DEM  production  process. 

Keywords:  Digital-Elevation-Model,  Geomorphometry,  Slope-Algorithm,  Aspect- 
Algorithm 


INTRODUCTION 

Digital  elevations  models  (DEMs),  rectangular  grids  of  elevation  values,  provide 
an  unparalleled  tool  for  geomorphometry.  DEMs  now  cover  a  wide  range  of 
scales,  from  global  through  continental  and  regional  to  local  and  almost  point 
scales.  Analysis  can  extract  a  wide  variety  of  parameters  relating  elevation,  slope, 
aspect  and  terrain  organization.  Conversely,  geomorphometry  provides  an  ideal 
tool  for  assessing  the  quality  of  DEMs.  Outliers  on  histograms  or  bivariate  plots 
of  geomorphometric  parameters  almost  always  represent  flawed  DEMs. 

At  small  scales  DEM  data  cover  the  entire  world,  but  at  large  scale  data 
availability  reflects  a  bias  toward  the  United  States.  This  bias  results  from  two 
factors:  (1)  the  United  States  government  has  one  of  the  most  aggressive 
programs  in  the  world  to  produce  DEMs,  and  (2)  almost  unique  in  the  world,  the 
United  States  national  mapping  agencies  place  no  restrictions  on  their  digital 
data.  DEMs  with  30  m  spacing  cover  almost  all  of  the  United  States,  and 
increasing  numbers  of  10  m  DEMs  are  coming  on  line.  With  data  spacing  ranging 
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from  almost  10  km  to  10  m,  available  DEMs  span  4  orders  of  magnitude. 
Experimental  1  m  resolution  LIDAR  (Light  Detection  and  Ranging)  DEMs  cover 
areas  almost  point-like  with  vast  quantities  of  data  (NOAA  Coastal  Services  Data 
Center,  2001),  and  will  add  another  order  of  magnitude  to  range  of  scales  when 
they  become  widely  available. 


DATA 

This  study  used  DEMs  from  continental  scale  (30"  spacing,  about  1  km), 
regional  scale  (3"  spacing,  60-90  m),  and  local  scale  (10-30  m  spacing).  Table 
1  summarizes  the  free  data  obtained  from  the  www.  These  data  sets  will  now  fit 
comfortably  on  hard  disks  of  mainstream  personal  computers.  All  of  these  data 
sets  use  binary  2  byte  elevations.  The  continental  data  sets  are  stored  uncompressed, 
whereas  the  other  data  sets  are  stored  with  the  same  standard  UNIX  compressed 
file  structure  obtained  from  the  USGS  web  site.  Because  compression  varies  with 
the  type  of  terrain  and  the  amount  of  smoothing  from  the  DEM  creation  process, 
the  estimates  below  should  be  viewed  only  as  rough  estimates  for  the  order  of 
magnitude  of  the  DEM  sizes.  Further,  if  future  large  scale  DEMs  use  4  byte  values 
for  elevation  as  the  USGS  is  starting  to  implement,  file  sizes  might  double. 

Three  different  United  States  government  mapping  agencies  have  created 
quasi-independent  continental  30"  DEMs:  DTED  Level  0  (NIMA,  2000),  GLOBE 
(NOAA  National  Geophysical  Data  Center,  2001),  and  GTOPO30  (USGS  EROS 
Data  Center,  2001).  After  downloading,  the  GLOBE  and  GTOPO30  datasets 
were  converted  to  1°  tiles  to  match  the  DTED,  and  provide  a  good  statistical 
sampling  unit  for  comparisons.  With  1°  elevation  tiles  each  data  set  will  fit  on  a 
single  CD-ROM,  if  some  of  the  interior  cells  in  Antarctica  are  removed.  The 
flatness  of  the  ice  sheet  and  the  convergence  of  the  meridians  make  the  30"  DEM 
a  poor  choice  for  the  polar  regions. 

The  250K  DEM  (USGS,  2001)  was  originally  created  by  a  precursor  of 
NIMA  and  contains  many  processing  anomalies  but  remains  a  valuable  regional 
data  resource  because  of  its  availability.  The  current  NIMA  product,  DTED  Level 
1,  provides  a  much  better  demonstration  of  what  this  data  resolution  can  achieve 
but  is  unfortunately  not  publicly  available.  The  DEMs  used  in  this  study  cover  all 
of  the  United  States  and  Puerto  Rico  except  for  Alaska.  The  SRTM  mission 


Table  1.  DEMs  available  for  comparison. 


Name 

DEMs  in  Sample 

Tile 

Spacing 

Sample 

Storage 

Full  US 

Full 

World 

ETOP05 

TerrainBase 

i 

5’ 

20  MB 

20  MB 

DTED0 

GLOBE, 

GTOPO30 

17,680 

21,803 

20,151 

1° 

30” 

650  MB 

650  MB 

USGS  250K 

940 

1° 

3” 

1  GB 

1  GB 

20  GB 

USGS  24K 

50,279  (42,921  unique) 

714’ 

30  m 

6.2  GB 

6.2  GB 

150  GB 

USGS  24K 

1066 

714’ 

10m 

1  GB 

50  GB 

1  TB 
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should  produce  this  world-wide  data  set  for  public  use  (Jet  Propulsion  Laboratory, 
2001),  and  the  first  samples  covering  some  of  the  United  States  appeared  in  early 
2002  and  South  America  in  the  summer  of  2003.  The  SRTMdata  appeared  too  late 
to  include  in  this  analysis,  and  will  require  careful  consideration  on  how  to  handle 
the  data  voids  when  conducting  statistical  analyses. 

Over  51,000  24K  DEMs  were  downloaded  from  the  same  USGS  site  along 
with  the  250K  DEM  (USGS,  2001)  in  the  summer  of  2000,  before  the  24K  DEMs 
were  removed.  USGS  allowed  an  automated  download  process,  the  only  practical 
way  to  obtain  such  a  large  sample.  This  data,  and  an  increasing  number  of  10  m 
DEMs,  must  now  be  downloaded  free  from  two  commercial  sites  that  restrict  the 
speed  or  number  of  DEMs  that  can  be,  and  which  forbid  the  use  of  automated  ftp, 
or  in  a  reprojected  format  from  the  National  Elevation  Dataset  (NED,  Gesch  and 
others,  2002).  While  the  data  considered  here  have  elevations  stored  as  2  byte 
integers  with  foot  or  meter  resolution,  some  newer  DEMs  now  use  2  byte  integer 
decimeters  and  others  use  4  byte  floating  point  values. 

The  24K  DEMs  have  two  characteristics  relevant  to  a  discussion  of 
geomorphometry.  The  “Level”  of  the  DEM  assesses  its  quality:  Level  1 ,  no  longer 
being  produced,  used  profiling  or  similar  photogrammetric  methods  to  produce 
lower  quality  DEMs,  whereas  the  newer  Level  2  DEMs  used  interpolation  from 
scanned  contour  lines  (USGS  National  Mapping  Division,  1997/1998).  Except 
for  contour  line  “ghosts”  where  elevations  corresponding  to  the  source  map 
contour  lines  are  overrepresented  in  the  DEM  (Guth,  1999c),  the  Level  2  DEMs 
provide  a  much  better  representation  of  the  terrain.  In  addition,  some  DEMs  have 
both  10  m  and  30  m  resolution.  Ten  meter  resolution  provides  much  greater  visual 
detail  in  representations  of  the  terrain. 

At  30  m  spacing  the  data  set  includes  42,921  (80%)  of  the  53,873  24K 
quadrangles  in  the  continental  United  States,  plus  81  DEMs  in  Puerto  Rico  and 
the  Virgin  Islands,  and  provides  an  excellent  sample  for  this  scale  data  at 
continental  resolution.  Hawaiian  data  were  not  available  on  the  USGS  web  site, 
and  Alaska  was  not  mapped  at  this  scale.  The  SRTM  mission  (Jet  Propulsion 
Laboratory,  2001)  should  produce  this  world-wide  data  set  for  US  military  use, 
although  the  radar’s  smoothing  might  make  the  effective  resolution  lower  than 
the  nominal  resolution.  The  data  set  includes  about  1000  10  m  level  2  DEMs,  or 
about  2%  of  the  total  area  of  the  continental  United  States. 

The  following  tests  on  the  differences  in  calculated  statistical  parameters 
will  be  conducted  using  these  DEMs  for  geomorphometry: 

1.  How  do  the  different  30"  continental  DEMs  differ? 

2.  How  do  the  30"  continental  and  3"  regional  DEMs  differ? 

3.  How  do  the  3"  regional  and  30  m  local  DEMs  differ? 

4.  How  do  the  Level  1  and  Level  2  local  DEMs  differ? 

5.  How  do  the  10  m  and  30  m  Level  2  local  DEMs  differ? 

GEOMORPHOMETRIC  PARAMETERS 

Pike  (1988)  introduced  the  concept  of  a  geometric  signature,  a  multi-variate 
description  of  topography  using  a  suite  of  measures,  and  later  expanded  the 
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concept  with  a  listing  of  49  variables  that  could  be  grouped  into  22  attributes 
(Pike,  2001).  He  considered  roughness  and  height  the  two  most  important 
attributes,  with  two  measures  of  texture  in  seventh  and  eleventh  position.  Fifteen 
different  variables  contribute  to  roughness. 

This  analysis  will  focus  on  three  geomorphic  parameters:  average  elevation, 
average  slope,  and  organization  strength.  These  parameters  reflect  the  interaction 
of  climate,  lithology,  and  topography.  Average  elevation,  combined  with  latitude, 
represents  a  primary  control  on  climate  for  both  precipitation  and  temperature. 
Slope  depicts  the  ruggedness  and  dissection  of  the  landscape.  Organization 
values  measure  the  processes  responsible  for  terrain  formation;  the  most  highly 
organized  regions  are  in  folded  mountain  belts  or  glacial  drumlin  fields.  Average 
elevation  and  average  slope  coincide  with  Pike’ s  (200 1 )  most  important  attributes, 
but  none  of  his  experimental  measures  fully  capture  this  organization  parameter. 

Computation  of  average  slope  requires  selection  of  a  slope  algorithm;  a 
number  have  been  proposed,  and  about  a  half  dozen  widely  used  (recent  summaries 
in  Skidinore,  1989;  Guth,  1995;  Hodgson,  1998;  Jones,  1998).  Results  for 
average  regional  slope  correlate  very  strongly  among  the  various  methods, 
although  the  absolute  values  vary.  This  analysis  uses  the  steepest  adjacent 
neighbor  algorithm,  but  comparisons  between  the  steepest  adjacent  neighbor  and 
four  nearest  neighbor  methods  show  how  the  choice  of  slope  method  affects 
geomorphic  results. 
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Fig.  1.  Average  slope  versus  flatness  parameter  ln(Sl/S2)  for  the  51,000  24K  USGS  DEMs.  The 
“good”  DEMs  show  the  high  correlation.  Outliers  consist  of  “bad”  DEMs  (visually  examined  to 
verify  processing  anomalies)  and  flat  and  overedge  DEMs  (identified  by  small  compressed  file 
sizes). 
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Eigenvector  methods  can  retrieve  a  flatness  parameter  and  organization 
parameter.  Chapman  (1952)  proposed  a  graphical  method  using  a  stereo  net  to 
characterize  terrain  orientations.  Woodcock  (1977)  described  an  eigenvalue 
method  for  fabric  shapes  in  structural  geology.  The  two  methods  can  be  combined 
for  automatic  characterization  of  terrain  organization  (Guth,  1999a,  1999b, 
2001).  Direction  consines  of  the  landscape  surface  were  computed  from  the  slope 
and  aspect  at  each  point,  and  then  eigenvectors  were  obtained.  The  eigenvalue 
method  extracts  two  independent  parameters  from  the  distribution  of  normals  to 
the  earth’ s  surface  calculated  from  the  DEM,  both  logs  of  the  ratios  of  eigenvalues. 
The  ratio  ln(Sl/S2)  reflects  a  flatness  parameter,  and  it  correlates  almost  perfectly 
with  the  average  slope  (Fig.  1).  The  correlation  is  negative  because  the  log  ratio 
is  flatness  and  not  steepness,  and  is  not  linear  because  of  the  log  ratio.  This 
correlation  should  not  be  surprising,  since  the  eigenvalue  method  starts  with  the 
slope  and  aspect  at  each  point  and  the  orientation  of  the  first  eigenvector  is  the 
average  normal  to  the  earth's  surface. 

The  ratio  ln(S2/S3)  reflects  the  strength  of  the  preferred  orientation  shown 
by  Fig.  2.  No  previously  defined  geomorphic  parameters  replicate  its  power  in 
quantifying  the  degree  of  organization  of  the  terrain.  Other  parameters  attempt  to 
measure  the  spatial  dimensions  of  periodic  landforms,  with  “grain”  representing 
the  distance  between  repeating  landforms.  Based  on  analysis  of  over  50,000  30 
m  DEMs,  Fig.  2  shows  the  three  most  highly  organized  quadrangles  in  the  United 
States.  The  figure  shows  the  dominant  direction  of  the  fabric,  defined  by  the 
orientation  of  one  eigenvector,  and  shows  how  it  swings  around  as  the  folded 
ridges  change  orientation  along  strike. 


Enders,  PA 
ln(S2/S3)  =  2.34 
orientation  =  69° 


Grantville,  PA 
ln(S2/S3)  =  2.30 
orientation  =  63° 


Clyman,  WI 


ln(S2/S3)  =  2.04 
orientation  =  187° 


Fig.  2.  The  three  most  highly  organized  quadrangles  in  the  United  States.  The  two  on  the  left  are 
adjacent  Level  1  DEMs  in  the  folded  Appalachians,  and  the  level  2  DEM  on  the  right  contains 
a  drumlin  field  in  Wisconsin.  One  km  UTM  grid  provides  scale. 
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Elevation  is  a  true  point  parameter,  defined  at  every  grid  location  in  the 
DEM.  While  slope  is  also  a  point  parameter,  it  requires  consideration  of  the 
neighboring  elevations  and  various  algorithms  use  from  2  to  8  adjacent  neighbors. 
The  most  commonly  used  algorithm  uses  four  neighbors  and  actually  ignores  the 
elevation  of  the  central  point  to  which  the  slope  is  assigned.  The  eigenvector 
method  uses  a  region  and  probably  requires  at  least  100  points  for  reasonable 
statistical  accuracy.  While  organization  can  be  calculated  as  a  point  parameter  for 
a  particular  region  size  about  a  point  with  a  sliding  window  over  the  DEM  (Guth, 
1999b,  2001),  this  study  will  concentrate  of  calculations  with  a  1°  region  size  for 
global  scale  or  7  1/2' region  for  continental  scale.  This  quadrangle  geomorphometry 
introduces  arbitrary  aspects  to  the  study,  but  sample  sizes  are  large  enough  that 
quadrangles  represent  reasonable  sampling  areas. 

METHODS 

All  analyses  used  the  MICRODEM  program  (Guth  and  others,  1987;  Guth, 
1995).  This  program  predates  commercial  GIS  programs  for  personal  computers, 
and  while  MICRODEM  has  now  become  a  full-feature  GIS  program  it  retains  a 
greater  emphasis  on  DEMs  and  geomorphometry  than  any  commercial  program. 
While  scripting  within  a  GIS  might  duplicate  the  operations  described  here  (e.g. 
contour  line  ghosts  or  terrain  organization),  I  prefer  embedding  the  capabilities 
within  MICRODEM  where  I  have  complete  control  of  the  operations  and 
displays.  With  object-oriented  programming,  the  DEM  object  contains  methods 
to  calculate  geomorphic  parameters  of  the  terrain. 

MICRODEM  indexes  the  DEMs,  and  sequentially  processes  them  to  calculate 
geomorphic  parameters.  For  the  30"  global  DEMs  the  program  used  1°  tiles.  For 
the  3"  250K  DEM,  the  comparison  used  1°  tiles  to  compare  with  the  30"  data,  and 
7  1/2'  tiles  to  compare  with  the  24K  DEM.  The  program  computed  five  variables 
for  each  tile;  average  slope,  average  elevation,  the  eigenvalue  flatness  and 
organization  parameters,  and  the  direction  of  terrain  organization.  The  database 
contains  five  additional  fields  for  the  DEM  name  and  its  bounding  geographic 
rectangle.  The  program  computed  28  additional  parameters  for  the  24K  DEM: 
duplicate  variables  to  compare  slope  algorithms,  DEM  quality  measures  like  the 
ghost  ratio  (Guth,  1999c),  and  additional  geomorphometric  parameters  which 
will  not  be  discussed  here  but  which  could  be  used  for  a  geometric  signature 
(Pike,  2001). 

The  GIS  database  files  for  the  geomorphic  parameters  allows  summary 
statistics,  bivariate  scatter  plots,  and  map  displays,  with  options  to  filter  the 
display — for  example,  to  only  show  terrain  organization  for  regions  with  greater 
than  a  specified  average  slope.  The  interactive  nature  of  the  GIS  allows  clicking 
on  a  point  on  the  graph  to  show  tabular  data  for  that  point  and  its  location  on  the 
map.  For  the  large  data  sets  considered  in  this  effort,  the  power  of  the  GIS 
represents  the  most  effective  way  to  explore  the  data. 

Drawing  on  earlier  work  of  Fenneman  dating  to  1917,  the  Fenneman  and 
Howard  (1946)  map  of  the  United  States  has  8  major  divisions,  25  provinces,  and 
86  sections  representing  distinctive  areas  having  common  topography,  rock  types 
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and  structure,  and  geologic  and  geomorphic  history.  A  digital  version  of  the  map 
(Hitt,  2002)  was  converted  to  geographic  coordinates,  and  then  each  quadrangle 
in  the  United  States  was  assigned  to  the  corresponding  section.  This  allows  the 
GIS  to  compute  statistics  based  on  the  Fenneman  physiographic  regions. 

Many  of  these  analyses  omit  24K  DEMs  identified  as  containing  “bad"  data, 
and  those  with  a  compressed  file  length  less  than  50  kb.  Small  file  sizes  occur  in 
two  cases:  coastal  “overedge”  DEMs  where  standard  quadrangles  leave  a  small 
bit  of  land  in  what  would  be  the  next  quadrangle,  and  extremely  flat  quadrangles 
in  regions  like  the  coastal  plains.  These  DEMs  account  for  a  disproportionate 
number  of  the  outliers  on  statistical  plots.  Overedge  DEMs  have  an  insufficient 
number  of  points  for  reliable  statistics,  and  the  DEMs  do  not  have  adequate 
vertical  resolution  for  good  statistics  in  very  flat  regions.  The  GIS  also  contains 
a  percent  ocean  field  for  the  global  data  sets;  often  the  coastal  DEMs  should  be 
ignored  for  the  eigenvector  analysis.  The  30"  1°  tiles  contain  only  about  1400 
elevations,  and  when  a  significant  proportion  are  missing  the  resulting  statistics 
represent  outliers  compared  to  full  data  sets. 

SLOPE  ALGORITHMS 

While  different  slope  algorithms  correlate  highly,  they  produce  different 
results.  Guth  (1995)  suggested  using  a  steepest  adjacent  neighbor  algorithm, 
while  the  most  popular  algorithm  appears  to  be  the  four  nearest  neighbors  (e.g. 
Hodgson,  1998;  Jones,  1998),  although  Evans  (1998)  makes  a  convincing  case  for 
the  superiority  of  an  eight  neighbors  algorithm.  Because  the  steepest  adjacent 
neighbor  and  four  nearest  neighbors  have  the  greatest  differences  among  the 
common  methods,  a  detailed  analysis  of  their  differences  highlights  the  effect  of 
slope  algorithm  on  geomorphic  parameters.  Because  slope  and  aspect  provide  the 
raw  input  to  the  eigenvector  algorithm,  the  slope  algorithm  might  affect  the 
organization  calculations. 

With  an  eight  point  neighborhood  about  the  central  point,  eight  partial  slopes 
can  be  calculated  in  each  of  the  eight  principal  cardinal  directions.  If  the  point 
elevations  are  considered  accurate,  these  represent  the  only  unambiguous  slopes 
that  can  be  determined.  The  challenge  for  the  slope  algorithm  is  to  use  these 
values  to  estimate  a  single  slope  at  the  central  point.  The  steepest  adjacent 
neighbor  algorithm  uses  the  steepest  of  the  eight  slopes:  this  slope  exists,  it  is 
measured  over  the  DEM  spacing  (or  the  diagonal),  and  in  many  applications  the 
steepest  slope  is  most  important.  The  four  nearest  neighbors  algorithm  uses  just 
four  neighbors:  it  is  measured  over  twice  the  DEM  spacing,  and  it  smooths  the 
results  by  using  twice  the  slope  distance  and  not  using  the  diagonal  neighbors. 

Figure  3  shows  a  comparison  of  the  average  quadrangle  slope  computed 
from  these  two  algorithms  for  30  m  DEMs  (28,479  level  2  and  18,689  level  1). 
This  analysis  omits  DEMs  identified  as  containing  bad  data,  and  those  with  a 
compressed  file  length  less  than  50  kb. 

Several  observations  are  clear  from  the  two  graphs:  (1)  the  two  slope 
algorithms  are  highly  correlated;  (2)  the  steepest  adjacent  neighbor  algorithm 
produces  consistently  larger  average  slopes;  (3)  the  level  2  DEMs  appear  to  be 
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Av-tfi  r*u*  Slops  ("i;,  sisepsst  adjacent  neighbor} 


Average  Slopa  (“&■,  iTnipnit  adjacent  neighbor) 


Fig.  3.  Average  slope  computed  by  two  different  algorithms,  for  both  Level  1  and  Level  2  USGS  24K 
DEMs.  DEMs  from  two  Fenneman  physiographic  provinces  indicated. 


steeper  than  the  level  1  DEMs,  with  both  algorithms  (this  is  not  proven  by  these 
graphs  because  the  DEMs  in  the  two  samples  are  not  the  same,  but  it  will  be 
proven  later);  (4)  the  wide  scatter  about  the  trend  lines  suggests  that  additional 
factors  affect  the  slopes. 
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Table  2.  Comparison  of  slope  algorithms  for  two  mountainous  quadrangles. 


Temple  of  Sinawava,  Utah 

Cascade  Pass,  WA 

Steepest 

Neighbor 

Four  Neighbors 

Steepest 

Neighbor 

Four  Neighbors 

Entire  DEM 

86.7% 

63.8% 

87.5% 

75.7% 

Peaks 

101.2 

33.4 

130.7 

35.8 

Other 

86.0 

64.6 

87.6 

75.7 

Table  3.  Comparison  of  correlation  coefficient  r  for  geomorphometric  parameters. 


DEM  scales  compared 

Average 

Elevation 

Average 

Slope 

Organization 

ln(S2/S3) 

DTEDO  vs  GLOBE,  n=17,559 

r=0.992 

r=0.975 

r=0.598 

DTEDO  vs  GTOPO30,  n=17,547 

r=0.992 

r=0.950 

r=0.516 

GLOBE  vs  GTOPO30,  n=20,l  14 

r=0.999 

r=0.979 

r=0.906 

DTEDO  vs  USGS  250K,  n=946 

r=1.000 

r=0.983 

r=0.680 

GTOPO30  vs  USGS  250K,  n=938 

r=0.999 

r=0.987 

r=0.701 

GLOBE  vs  USGS  250K,  n=941 

r=0.999 

r=0.987 

r=0.670 

USGS  250K  vs  24K,  n=5 1 ,345 

r=0.997 

r=0.951 

r=0.496 

USGS  24K,  Level  2,  10  vs  30  m,  n=472 

r=1.000 

r=0.998 

r=  0.960 

USGS  24K,  Level  1  vs  2,  30  m,  n=6476 

r=  1.000 

r=0.968 

r=0.545 

To  investigate  the  effect  of  additional  factors  on  the  computed  average 
slopes,  two  level  2  DEMs  were  selected  from  the  steeper  side  of  the  graph.  Both 
had  similar  average  slopes  using  the  steepest  adjacent  neighbor  method,  but 
greatly  different  slopes  for  the  four  adjacent  neighbors  method  (Table  2). 

As  noted  by  Guth  (1995),  these  two  algorithms  will  differ  most  for  peaks, 
pits,  ridges  and  valleys.  Pits  are  extremely  rare,  both  in  nature  and  the  DEMs,  but 
the  other  terrain  types  occur  about  in  5-6%  of  the  points  in  these  DEMs.  While 
the  differences  in  these  points  can  be  huge  (see  the  results  for  peaks  in  Table  2), 
it  is  the  other  points  that  make  the  difference — the  95-96%  of  the  points  classified 
as  “other”  have  average  slopes  almost  identical  to  the  averages  for  the  entire  DEM 
(Table  3).  On  this  measure  these  two  DEMs  differ  significantly. 

Analysis  of  the  graphs  in  Fig.  3  revealed  that  the  high  slope  DEMs  occurred 
primarily  in  two  geomorphic  provinces,  the  Grand  Canyon  and  the  Northern 
Cascades.  A  few  equally  steep  DEMs  occur  in  other  provinces,  including  the 
Temple  of  Sinawava  DEM  in  the  High  Plateaus  of  Utah  Section  adjacent  to  the 
Grand  Canyon.  Highlighting  those  provinces  on  the  graph  showed  two  very 
distinct  trends:  the  four  adjacent  neighbors  method  produces  gentler  slopes  in  the 
Grand  Canyon  than  in  the  Northern  Cascades,  for  DEMs  that  the  steepest  adjacent 
neighbor  method  considered  equally  steep.  Figure  4  shows  shaded  relief  depictions 
of  these  two  DEMs.  The  Temple  of  Sinawava  quadrangle  has  more  youthful 
dissection.  More  points  have  extreme  slopes  along  one  of  the  diagonals,  and  these 
influence  the  average  slope  of  the  quadrangle  only  with  the  steepest  adjacent 
neighbor.  The  commonly  used  four  neighbors  slope  algorithm  behaves  differently 
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Fig.  4.  Shaded  relief  views  to  Level  2,  30  m  USGS  DEMs  of  Temple  of  Sinawava  and  Cascade  Pass. 


Organization  Parameter  (Steepest  neighbor) 


Fig.  5.  Organization  parameter  for  30  m  Level  2  USGS  DEM  calculated  with  two  different  slope 
methods. 


HANGING  ROCK  CANYON.  CA  •  24000  Slope  0-52  HANGING  ROCK  CANYON.  CA  •  24000  Slope  5-102  HANGING  ROCK  CANYON.  CA  •  24000  Slope  10-202 

r>*8204  n=13477  n»27665 
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Fig.  6.  Aspect  distributions  for  the  Level  2  Flanging  Rock  Canyon  30  m  DEM,  divided  by  slope  class. 
This  DEM  has  elevations  in  meters. 


TROPIC  CANYON.  UT-24000  Slope  0-5X  TROPIC  CANYON.  UT -24000  Slope  5-10*  TROPIC  CANYON.  UT-24000  Slope  10-20* 
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DEM  has  elevations  in  decimeters. 
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in  different  physiographic  provinces,  compared  to  the  steepest  adjacent  neighbor 
algorithm. 

Since  the  eigenvector  technique  relies  on  the  slope  as  input,  Fig.  5  shows  the 
organization  parameter  for  level  2  DEMs  calculated  with  the  two  slope  algorithms. 
Scatter  is  minimal,  the  graph’s  slope  is  close  to  1,  and  the  differences  do  not 
correlate  with  the  physiographic  province.  The  eigenvector  technique  appears  to 
extract  organization  independent  of  the  slope  algorithm  use. 

ASPECT  ALGORITHM 

Selection  of  an  aspect  algorithm  presents  much  less  choice  than  that  of  a 
slope  algorithm,  because  extending  many  of  the  slope  algorithms  would  result  in 
a  very  limited  number  of  possible  aspect  directions.  Guth  (1995)  suggested  use 
of  the  eight  nearest  neighbors  with  even  weighting,  and  that  algorithm  has  been 
used  here.  However,  this  algorithm  does  not  provide  an  even  or  uniform  distribution 
of  aspects. 

Figure  6  shows  the  aspects  for  a  30  m  DEM.  Note  the  spikes  in  the  aspect 
distribution,  especially  prominent  for  the  low  slope  parts  of  the  DEM.  Aspects  at 
45°  intervals  are  overrepresented,  and  the  algorithm  appears  to  work  well  only 
when  the  slope  exceeds  30%.  The  problem  originates  in  the  discrete  nature  of  the 
elevations:  with  meter  elevation  resolution  and  30  m  spacing,  the  basic  slope 
values  only  change  by  1/30.  To  verify  this  assumption.  Fig.  7  shows  the  aspect 
distributions  for  a  10  m  DEM  with  elevations  in  decimeters.  This  allows  slope 
values  to  change  by  1/100,  and  in  this  case  the  algorithm  appears  to  produce 
uniform  and  natural  aspect  distributions  if  the  slope  is  steeper  than  10%. 

RESULTS 

Table  3  summarizes  the  results  of  the  comparisons  of  DEMs  at  different 
scales  or  sources.  The  correlation  coefficients  for  both  average  slope  and  average 
elevation  exceed  0.95  for  all  the  comparisons,  but  with  these  sample  sizes  such 
correlation  coefficients  still  allow  some  scatter  and  some  significant  anomalies. 
The  organization  parameter  shows  more  scatter  and  lower  correlations. 

1.  Continental  data  sets 

The  continental  data  sets  agree  with  each  other  extremely  well  for  average 
slope  and  average  elevation,  and  less  well  for  organization  level.  Table  3  shows 
correlation  coefficients  greater  than  0.99  for  elevation  and  0.95  for  slope,  and  as 
low  as  0.51  for  organization.  GTOPO30  and  GLOBE  show  the  strongest 
correlation,  and  both  have  complete  coverage  of  the  world.  DTED  level  0  does  not 
include  Antarctica,  and  still  has  a  large  gap  in  the  Amazon  Basin. 

Despite  correlation  coefficients  for  average  elevation  exceeding  0.99, 
differences  in  average  elevation  for  1°  cells  can  exceed  2000  m.  Six  cells  in  the 
Andes  of  Peru  and  Bolivia  have  differences  this  large  between  DTED  and 
GLOBE,  and  one  has  a  difference  exceeding  3000  m.  This  suggests  that  for  parts 
of  the  world  the  continental  data  sets  still  have  significant  problems. 
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Fig.  8.  Comparison  of  the  organization  parameter  ln(S2/S3)  calculated  from  DTED  Level  0  and 
GLOBE  for  1°  cells.  Note  that  (a)  has  a  different  scale  than  the  others. 
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To  investigate  how  well  the  continental  data  sets  work  in  regions  with  good 
digital  representations  of  topography.  Fig.  8  shows  the  scatter  plots  for  organization 
plotted  from  DTED  and  GLOBE.  Figure  8a  shows  the  global  data,  with  a  great 
deal  of  scatter  and  a  correlation  of  0.594.  With  the  data  restricted  to  part  of  North 
America  (N23-N54,  and  W130-W73),  Fig.  8b  shows  much  less  scatter  and  an 
improved  correlation  of  0.887.  Restricting  the  data  to  regions  with  slopes  greater 
than  2%  and  5%  in  Figs.  8c  and  8d  increases  the  correlation  even  further  to  0.97. 
Organization  values  have  the  least  reliability  in  regions  of  low  slope. 

Assessing  the  accuracy  of  DTED  Level  0  versus  GLOBE  or  GTOPO30 
remains  beyond  the  scope  of  this  work,  because  it  likely  varies  regionally.  The 
tremendous  differences  in  average  elevation  noted  above  occur  on  the  eastern 
flank  of  the  Andes,  where  DTED  shows  the  average  elevation  in  one  cell  to  be 
over  2000  m  higher  than  the  cell  immediately  to  the  east,  and  neighboring  cells 
do  not  even  remotely  match  along  their  common  edge.  Only  a  few  of  the  cells  in 
that  region  used  higher  resolution  DTED  as  their  source,  and  most  of  the  others 
show  only  the  grossest  outlines  of  topography.  GTOPO30  and  GLOBE  improve 
the  elevation  picture  somewhat  in  that  region,  but  clearly  show  which  cells  had 
high  resolution  source  data  available.  The  SRTM  30"  global  set  has  just  been 
released  (summer  2003),  and  should  now  be  used  for  this  scale  analysis. 

2.  Continental  versus  regional  DEMs 

This  comparison  must  be  restricted  to  the  United  States  because  of  the  lack 
of  publicly  available  data  elsewhere.  Table  3  shows  the  comparison  of  about  940 
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Fig.  9.  Comparison  of  the  organization  parameter  between  DTED  0  and  the  USGS  3",  250K  DEM. 
This  shows  the  432  cells  covering  the  United  States  with  an  average  slope  in  the  USGS  DEM 
greater  than  5%. 
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Table  4.  USGS  250K  versus  24K,  7  1/2'  tiles. 


Average  slope 

Number  ofl'A'  tiles 

Organization  Parameter 
Correlation  Coefficient  r 

All  data 

34,887 

0.473 

>2% 

23,241 

0.612 

>5% 

13,026 

0.757 

>10% 

7,956 

0.818 

>20% 

3,598 

0.818 

>30% 

1496 

0.835 

1°  tiles  for  DTED  0,  GLOBE,  and  GTOPO30  with  the  USGS  250K  DEM.  The 
three  continental  data  sets  have  slightly  different  numbers  of  tiles  because  of 
differences  in  continental  masking.  Average  elevation  and  average  slope  correlate 
almost  perfectly,  showing  that  the  continental  DEMs  accurately  capture  this 
aspect  of  topography. 

Organization  correlates  less  well,  with  correlation  coefficients  around  0.68. 
As  with  the  comparison  among  the  continental  DEMs,  the  disagreements  tend  in 
occur  in  the  flatter  DEMs.  For  the  DTED  0  comparison  with  the  USGS  3"  250K 
DEM,  the  correlation  coefficient  increases  from  0.6802  (n  =  946)  for  all  data,  to 
0.7314  for  DEMs  with  an  average  slope  greater  than  2%  (n  =  642),  and  to  0.7875 
for  DEMs  with  an  average  slope  greater  than  5%  (n  =  432,  Fig.  9).  Areas  of 
different  slope  will  occur  within  1°  cells,  so  that  low  slope  areas  will  produce 
different  estimates  of  organization  when  aggregated  over  the  entire  cell. 

5.  Regional  versus  local  DEMs 

Table  3  shows  that  the  USGS  250K  and  24K  DEMs  show  almost  perfect 
correlations  in  average  elevation  and  slopes.  At  the  quadrangle  level,  despite 
severe  contour  line  ghosts,  the  250K  DEM  accurately  depicts  average  terrain 
characteristics.  The  organization  parameter  between  the  two  series  correlates  less 
well,  but  as  Table  4  shows,  the  flat  DEMs  account  for  most  of  the  poor 
correlations. 

4.  Scale  in  local  DEMs  10  m  versus  30  m 

Somewhat  surprisingly.  Tables  3  show  that  for  these  geomorphic  parameters 
the  10  m  DEMs  do  not  provide  increased  precision.  The  10  m  DEMs  are  visually 
vastly  superior  in  the  fine  details  they  provide,  but  produce  essentially  identical 
statistical  results.  Even  the  organization  parameter  shows  a  high  correlation 
coefficient,  and  it  only  increases  from  0.954  to  0.981  as  average  slope  increases 
from  2%  to  greater  than  30%. 

5.  Quality  in  local  DEMs 

USGS  produces  two  levels  of  DEMs,  and  in  this  sample  they  were  present  in 
almost  equal  numbers.  Table  5  shows  how  the  calculated  organization  parameters 
correlate  as  the  average  slope  in  the  DEM  increases. 
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Table  5.  USGS  24K,  7  1/2'  tiles,  Level  1  versus  Level  2. 


Average  slope 

Number  of  IVi  tiles 

Organization  Parameter 
Correlation  Coefficient  r 

All  data 

6476 

0.545 

>2% 

6454 

0.544 

>5% 

5940 

0.606 

>10% 

4695 

0.706 

>20% 

2909 

0.757 

>30% 

1586 

0.890 
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Fig.  10.  Shaded  reflectance  maps  of  Level  1  and  Level  2  DEMs  of  the  same  quadrangle  in  New 
Mexico. 


One  reason  for  the  scatter  in  this  comparison  is  the  wide  range  in  quality  of 
the  level  1  DEMs.  Figure  10  shows  the  level  1  and  level  2  DEMs  for  a  very  flat 
region,  and  the  Level  1  DEM  is  one  of  the  worst  examples  available.  Its  calculated 
organization  parameter  would  rank  it  as  one  of  the  most  organized  in  the  country, 
but  as  the  Level  2  DEM  shows,  the  organization  is  entirely  an  artefact  of  the  DEM 
creation  process. 

6.  Average  slope  at  different  DEM  scales 

As  has  been  determined  in  previous  studies,  computed  slope  increases  as  the 
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Table  6.  Comparison  of  average  slope  from  different  scale  DEMs. 


Small  Scale  DEM  and  Tile  Size 

Sample  Size 

Slope  Ratio  to  Large  Scale  DEM 

GLOBE  Sloped0  tiles) 

n=941 

78%  250K  slope 

250K  Slope  (TA’  tiles) 

n=5 1,345 

86%  24K  Slope 

24K  level  1  Slope  (714’  tiles) 

n=6476 

85%  24K  Level  2  slope 

24K  (30  m)  level  2  Slope  (714’  tiles) 

n=472 

100%  24K  (10  m)  level  2  Slope 

Average  slope  level  1  DEM 

Fig.  1 1 .  Comparison  of  the  average  slope  (steepest  adjacent  neighbor  method)  for  6226  30  m  DEMs 
with  both  Level  1  and  level  2  coverage. 


size  of  the  sampling  region  decreases.  This  study  quantifies  that  relationship, 
using  large  samples  of  real  DEMs.  Table  6  shows  the  relationship  observed,  and 
Fig.  1 1  shows  the  relationship  between  independent  level  1  and  level  2  DEMs  of 
the  same  quadrangle.  Much  of  the  scatter  is  likely  due  to  a  minority  of  very  poor 
quality  Level  1  DEMs  (such  as  that  in  Fig.  10).  The  only  surprise  was  that  10  m 
and  30  m  DEMs  produced  identical  average  slopes. 

7.  DEM  quality  control 

First  attempts  to  create  bivariate  plots  of  geomorphic  parameters  from  the 
USGS  24K  DEMs  led  to  the  recognition  that  (1)  overedge  and  extremely  flat 
DEMs  produced  Door  statistics,  and  (2)  most  outliers  on  the  plots  represent 
problems  in  the  source  DEMs.  This  problem  was  severe  enough  that  a  “B  AD_DEM” 
field  was  added  to  the  database,  and  a  quick  way  to  flag  these  DEMs  was 
incorporated  into  MICRODEM.  Figure  10  shows  one  particularly  bad  DEM, 
identified  with  striping  (Oimoen,  2000)  when  it  was  visually  inspected  as  a 
candidate  for  the  most  organized  DEM  in  the  United  States. 
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Table  7.  “Bad”  24K  DEMs  by  series. 


DEM  Series 

Number  Available 

Number  “Bad” 

Percent  “Bad” 

Level  1,  30  m 

19,444 

256 

1.31% 

Level  2,  30  m 

30,835 

68 

0.22% 

Level  2,  10  m 

1066 

6 

0.56% 

Average  quadrangle  slope  initially  revealed  the  most  defective  DEMs, 
because  the  bad  elevations  are  usually  so  much  larger  than  the  others  that  they 
affect  the  average.  Figure  1  shows  the  strong  relationship  between  average 
quadrangle  slope  and  the  flatness  parameter.  After  first  flagging  flat  and  overedge 
DEMs,  a  number  of  outliers  remained.  These  had  a  higher  average  slope  than  their 
flatness  parameter  would  indicate.  Clicking  the  point  on  the  graph  brings  up  the 
database  table  for  that  record,  and  a  click  there  opens  the  DEM.  When  visual 
analysis  confirms  that  the  DEM  is  defective,  it  can  be  flagged  and  ignored  for 
further  analysis.  This  pair  of  variables  works  because  bad  values  affect  the 
arithmetic  mean  used  for  average  slope  more  than  they  affect  the  vector  average 
used  for  the  flatness  parameter. 

The  following  errors  have  been  identified:  (1)  large  rectangular  regions  of 
the  DEM  with  uniform,  erroneous  elevations  (probably  a  flag  representing 
regions  of  no  data,  appropriate  to  processing  software  but  not  to  the  SDTS  DEM 
profile);  (2)  a  single  partial  row  or  column  of  bad  elevations  along  one  margin  of 
the  DEM  (often  one  or  two  bad  postings);  (3)  a  single  spike  elevation,  in  some 
cases  where  a  posting  in  feet  was  inserted  into  a  DEM  in  meters;  (4)  snaking 
“lines”  of  bad  elevations  within  the  DEM;  (5)  DEMs  mislabeled  as  feet  or  meters 
instead  of  the  correct  value,  which  show  up  in  regional  maps  of  average  elevation 
where  they  are  either  three  times  or  one  third  the  elevation  of  their  neighbors;  (6) 
lakes  with  a  uniform  elevation  that  does  not  match  the  surroundings,  or  the  height 
on  the  topographic  map.  The  striping  artefacts  in  some  level  1  DEMs  (Oimoen, 
2000)  are  not  included,  because  I  have  not  been  able  to  automate  their  identification, 
and  assessing  their  severity  remains  a  very  subjective  process. 

At  present  at  least  330  of  the  5 1 ,345  DEMs  downloaded  from  the  USGS  have 
been  flagged  as  defective;  this  is  a  minimum  estimate  because  only  the  most 
obvious  statistical  outliers  have  been  examined  visually  to  confirm  the  DEM 
quality.  After  redownloading  a  sample  of  these,  it  is  clear  the  files  posted  on  the 
web  were  defective  and  not  the  download  process.  With  compressed  downloads, 
a  faulty  download  is  more  likely  to  produce  a  file  that  will  not  decompress  rather 
than  one  with  erroneous  elevations.  Table  7  suggests  that  over  1%  of  the  Level 
1  DEMs  have  errors,  and  that  for  Level  2  DEMs  the  error  rate  has  decreased.  This 
analysis  may  not  reflect  the  current  status  of  the  USGS  data  because  (1)  the  USGS 
performed  a  major  reprocessing  effort  in  the  summer  of 2001  to  correct  registration 
problems  with  some  DEMs,  (2)  increasing  numbers  of  10  m  DEMs  are  now 
available,  and  (3)  it  is  not  currently  feasible  to  duplicate  collection  of  the  50,000 
DEM  data  set. 
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Maximum  point  slope  is  probably  the  single  best  parameter  for  flagging  bad 
DEMs.  While  a  few  DEMs  have  valid  huge  point  slopes,  like  the  volcanic  neck 
at  Ship  Rock,  New  Mexico  (slope  1561%)  or  the  cliff  face  of  Half  Dome  in 
Yosemite  National  Park  (slope  1553%),  excessive  point  slope  values  over  1000% 
almost  always  represent  bad  postings  in  the  DEM  (at  least  29 1  of  364  bad  DEMs 
in  this  sample).  While  many  of  the  errors  in  these  DEMs  are  not  severe  and  would 
not  significantly  affect  most  uses  of  the  DEM,  they  could  be  eliminated  with  a 
simple  quality  control  step.  The  DEM  should  be  displayed  as  a  shaded  reflectance 
map,  the  location  of  the  steepest  point  slope  flagged,  and  the  elevations  of  the 
postings  surrounding  the  point  displayed.  The  maximum  slope  could  be  compared 
to  the  distribution  of  maximum  slope  with  the  Fenneman  province,  and  if  the 
maximum  slope  in  the  DEM  occurs  in  a  pit,  the  point  almost  certainly  represents 
an  error.  The  DEM  could  then  be  accepted,  or  additional  edits  performed.  While 
this  sample  was  downloaded  in  the  summer  of  2000,  some  of  the  bad  DEMs  are 
still  present  on  the  commercial  servers  now  distributing  the  USGS  data.  Osborn 
and  others  (2001)  describe  the  challenges  faced  by  USGS  in  developing  quality 
control  methods  for  the  DEMs. 


CONCLUSIONS 

•  Different  slope  algorithms  behave  differently  depending  on  the 
physiographic  province.  Most  of  the  difference  between  the  steepest  adjacent 
neighbor  algorithm  and  the  four  nearest  neighbors  algorithms  is  not  due  to  the 
small  number  of  peaks,  ridges,  and  valleys,  which  can  have  huge  differences,  but 
to  the  vast  majority  of  points  which  cannot  easily  be  categorized  over  a  nine  point 
region. 

•  Aspect  algorithms  do  not  work  well  in  low  relief  regions,  producing  too 
many  aspects  at  45°  intervals.  They  never  work  well  at  very  low  slopes,  and  the 
steepness  at  which  they  begin  to  produce  uniform  results  appears  to  be  related  to 
the  ratio  of  vertical  resolution  to  horizontal  resolution  in  the  DEM. 

•  Different  scales  of  DEMs  accurately  capture  average  elevation,  which 
correlates  very  strongly  across  scales. 

•  Different  scales  of  DEMs  accurately  capture  average  slope,  which 
correlates  very  strongly  across  scales.  The  absolute  value  of  average  slope 
increases  as  the  resolution  of  the  DEM  increases  (post  spacing  decreases), 
although  there  appears  to  be  no  change  between  data  spacings  of  10  and  30  m. 

•  Different  scales  of  DEMs  accurately  capture  terrain  organization,  both 
in  orientation  and  magnitude.  The  correlation  is  weakest  in  low  slope  regions, 
probably  because  of  problems  with  the  aspect  algorithms,  and  is  strongest  in 
moderate  to  high  relief  areas. 

•  Geomorphic  parameters  can  provide  quality  control  for  DEMs. 
Experience  with  the  USGS  DEMs  suggests  that  a  simple  maximum  point  slope 
check  should  be  performed  on  all  DEMs,  with  excessive  values  over  1000% 
almost  always  representing  bad  postings  in  the  DEM.  This  procedure  should  be 
implemented  in  the  DEM  production  process. 
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